Electron Correlation via Frozen Gaussian Dynamics 
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We investigate the accuracy and efficiency of the semiclassical Frozen Gaussian method in describ- 
ing electron dynamics in real time. Model systems of two soft-Coulomb-interacting electrons are used 
to study correlated dynamics under non-perturbative electric fields, as well as the excitation spec- 
trum. The results show that a recently proposed method that combines exact-exchange with semi- 
classical correlation to propagate the one-body density-matrix holds promise for electron dynamics 
in many situations that either wavefunction or density-fimctional methods have difficulty describing. 
The results also however point out challenges in such a method that need to be addressed before it 
can become widely applicable. 



I. INTRODUCTION 

Accurately capturing the correlated motion of elec- 
trons in atoms, molecules, and solids remains an active 
research area today. In ground-state electronic structure 
problems, the difficulties of solving many-electron sys- 
tems are rolled up into the correlation energy, a.k.a. the 
stupidity energy [1]. It has been examined in many dif- 
ferent scenarios and various limits, but its complexity 
still haunts us today. Essentially representing the devi- 
ation of the many-electron wavefunction from an anti- 
symmetrized product of single-particle orbitals, corre- 
lation plays a significant role in time-dependent prob- 
lems as well. When atoms and molecules are exposed to 
either perturbative or non-perturbative external fields, 
fascinating and subtle electron interaction effects mean 
one must go beyond the "single active electron" pic- 
ture, yet to solve the time-dependent Schrodinger equa- 
tion (TDSE) for more than two electrons in strong fields 
pushes today's computational limits [2, 3]. 

The advent of time-dependent density functional the- 
ory (TDDFT) in 1984 [4] opened up the possibility of 
using single-particle orbitals to describe the dynam- 
ics of interacting electrons exactly. In TDDFT a sin- 
gle slater determinant (SSD) in a non-interacting sys- 
tem (the Kohn-Sham (KS) system) is propagated in time, 
such that it reproduces the exact time-dependent den- 
sity of the true interacting system. Clearly the KS deter- 
minant is not the true interacting wavefunction, nor is it 
supposed to be an approximation to it; nevertheless by 
the rigorous theorems of TDDFT all observables of the 
true, correlated, electronic system can be in principle ex- 
tracted from it exactly In practise, approximations are 
needed for the exchange-correlation functionals, limit- 
ing the accuracy of the method. TDDFT is a method of 
choice in the linear-response regime where the photo- 
spectra of atoms, molecules and clusters can be found 
accurately The computational effort involved scales rel- 
atively well, so the spectrum of large systems such as 
the green fluorescent protein[5], or even candidates for 
solar cells [6], have been studied. 

Although it has proved successful for a wide vari- 
ety of time-dependent problems [7, 8], its progress for 



real-time dynamics in non-perturbative fields has been 
slower. Three major obstacles involve the lack of mem- 
ory in the commonly-used functional approximations, 
lack of good functional approximations to extract ob- 
servables not directly related to the density, and the in- 
accuracy of the usual functional approximations when 
the true wavefunction evolves far from a SSD [9, 10]. 

Recently a method that accounts for electron corre- 
lation semiclassically has been proposed [10] that will, 
in principle, remedy the aforementioned problems of 
TDDFT. The method operates within the context of time- 
dependent density-matrix functional theory (TDDMFT), 
where the one-body density matrix is propagated in 
time. The second-order density-matrix enters into the 
equation for the one-body density-matrix, so is approxi- 
mated as a functional of the one-body density-matrix in 
TDDMFT. However it is still difficult to find approxima- 
tions in TDDFMT that work well when the true wave- 
function evolves from being close to a SSD to being 
far from one [11-13], i.e. those developed so far can- 
not dynamically change "occupation numbers". In the 
proposed method of Ref. [10], a separate semiclassical 
propagation is done from which the correlation poten- 
tial is extracted, and used to drive the actual propaga- 
tion. The approach naturally contains memory effects 
carried along by classical trajectories, all one-body ob- 
servables may be immediately extracted, and occupa- 
tion numbers do dynamically evolve. 

Whether the proposed method is accurate for a given 
application however depends on how close the semi- 
classical dynamics is to the true dynamics. In the 
method, only the correlation component of the dy- 
namics is treated semiclassically, while all other terms 
(including Hartree and exchange), are treated exactly. 
However, if the semiclassical dynamics deviates far 
from the true d5mamics, the proposed method is un- 
likely to be accurate. Therefore it is worth investigat- 
ing how well a semiclassical description of the entire dy- 
namics is for the problems of interest. This is the goal 
of the present paper. We compute the semiclassical dy- 
namics (within the Frozen Gaussian approximation) of 
model two-electron systems where the exact dynamics 
can be numerically exactly computed, so that we can 



test the accuracy of the semiclassical dynamics against 
this. We can also compare with the exact KS system: 
that is, using the exact exchange-correlation potential 
(that yields the exact density), but comparing observ- 
ables such as momentum-densities, computed directly 
from the KS SSD (instead of using the unknown appro- 
priate observable-functionals as in exact TDDFT). Our 
examples model the following phenomena: states of 
double-excitation in spectra, dynamics in a strong os- 
cillating field, and population transfer to excited states 
via resonant driving, or via an optimal field. The chosen 
systems are one-dimensional models of a quantum dot 
and a helium atom. 

We show both successes and failures of the semiclassi- 
cal propagation for these applications. We expect how- 
ever that the proposed method of Ref. [10] is more ac- 
curate than such a semiclassical treatment of the entire 
problem since there it is only the correlation component 
that is approximated semiclassically. So the results of 
our studies are expected to improve once applied in the 
context of the method of Ref. [10]. 

Our results also have relevance for semiclassical stud- 
ies of djmamics in their own right (i.e. outside of 
density-functional- type methods). Although semiclassi- 
cal treatments of dynamics have seen much interest for 
nuclear dynamics (eg. Refs. [14-16]), there have only 
been a few applications to electron dynamics [17-23] 
and almost all involved a single electron. Part of the rea- 
son is that in atoms there is a significant probability that 
classical trajectories will autoionize after a few cycles in 
the atomic well: through the electron interaction, one 
electron gains energy from the other and ionizes while 
the other drops below its zero-point energy; a process 
that is classically allowed but quantum-mechanically 
forbidden. This is a manifestation of what has been 
called the "zero-point-energy" problem [23-25], and 
happens also when there are coupled vibrational modes, 
for example. These trajectories create a lot of noise in 
the semiclassical sum, making the convergence very dif- 
ficult, and so are typically terminated. For example, 
in Ref. [23], accurate spectra of the He atom in both 
its collinear Zee and eZe configurations were obtained, 
by discarding trajectories where one electron reaches a 
certain threshold distance. Theoretically, a question is 
whether fundamentally these trajectories should be ex- 
cluded from the semiclassical sum. That is, if the semi- 
classics could be done with an infinite number of trajec- 
tories, whether these trajectories incorrectly contribute 
to the semiclassical sum. We find that indeed their con- 
tribution to the semiclassical sum decreases as more tra- 
jectories are added: in this sense, they are valid classical 
trajectories, and via phase-interference, the semiclassi- 
cal dynamics restores the zero-point-energy, eliminating 
their contribution from the semiclassical sum. In prac- 
tise, nevertheless, we must deal with a finite number of 
trajectories, so sensible ways to discard these trajectories 
must be devised. Our results on the model atoms, as 
well as those on a model quantum dot, lend support to 



the conclusions of Ref. [23] that semiclassical dynamics 
is useful for electrons, for both spectra and dynamics in 
non-perturbative fields (in some cases). A draw of semi- 
classics is the interpretive power it carries and we hope 
to use it also to to interpret mechanisms for processes in- 
volving interacting electrons, e.g. the role of correlation 
in multiphoton ionization and multielectron ionization. 

The "flavor" of semiclassics explored in the current 
work is the Frozen Gaussian (FG) dynamics, originally 
proposed by Heller [26]. This is not one of the rigorous 
semiclassical propagators that satisfy the TDSE up to or- 
der fi; it is however closely related to the rigorous Heller- 
Herman-Kluk-Kay (HHKK) propagator [14-16, 29-31], 
which is based on a coherent-state representation of the 
semiclassical propagator. In Section II, after a brief re- 
view of TDDFT and TDDMFT that explains the moti- 
vation behind Ref. [10], we review some background 
regarding FG dynamics and how it is proposed to be 
used in the method of Ref. [10] in density-matrix prop- 
agation. Section III presents results on one-dimensional 
model systems of two electrons, using FG for dynam- 
ics and spectra. Finally, we make some conclusions in 
Sec. IV. 



II. BACKGROUND 

In this section, we briefly review TDDFT, the moti- 
vation behind going from TDDFT to a density-matrix 
approach, why we consider a semiclassical approach to 
correlation, and semiclassical dynamics. 



A. TDDFT and its challenges 

TDDFT is an exact reformulation of the non- 
relativistic time-dependent quantum mechanics of 
many-body systems [4, 8]. It operates by mapping the 
true problem of interacting electrons into a ficitious, 
non-interacting system of fermions, called the Kohn- 
Sham (KS) system, whose one-body density is, in prin- 
ciple, exactly that of the true system, and from which, 
in principle, all properties of the true system can be ob- 
tained. The many-body effects of the interacting sys- 
tem are modelled by a one-body potential, called the 
exchange-correlation potential Uxc- Similar thus in fla- 
vor to ground-state density functional theory, its func- 
tional are nevertheless quite different. For example, Uxc 
functionally depends on the entire history of the den- 
sity, as well as the initial-state of the true system and 
the initial state of the KS system, v^c {n, *o, *o] (r, t)[32]. 
This memory-dependence is however lacking in most of 
the applications today: they use an "adiabatic" approx- 
imation, which simply feeds the time-evolving density 
into a ground-state approximation: 

<^-[n,vI/o,$o](rt) = <cW(r)U(r)=n(rt) . (1) 
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With this simple approximation, TDDFT has had great 
success in calculating spectra and response: com- 
putationally it scales similar to methods such as 
time-dependent Hartree-Fock (TDHF) or configuration- 
interaction singles, but with accuracy that is usually 
far greater. Yet, this simple approximation is also 
behind why in some cases approximate TDDFT fails: 
e.g. to capture states of double-excitation charac- 
ter a frequency-dependent kernel is required [33-35], 
but the adiabatic approximation yields a frequency- 
independent kernel. 

The lack of memory in the usual functional approx- 
imations is also one of the reasons why the applica- 
tion of TDDFT to real-time dynamics applications has 
not progressed so fast. Model systems for which exact 
results are available indicate that memory-dependence 
can sometimes be a crucial factor in the exchange- 
correlation potential [36-40]. Another difficulty in 
real-time dynamics, is the problem of "observable- 
functionals": although in principle all properties of in- 
terest can be extracted from the KS system, it is not 
known hozu to extract those that are not directly re- 
lated to the density. One may write these observables 
as operators on the exact wavefunction, however sub- 
stituting the KS wavefunction is usually a poor ap- 
proximation. For example, for momentum distribu- 
tions, it was shown[9] for a model system with one 
electron ionizing that the KS momentum distribution 
incorrectly develops strong oscillations as the electron 
moves away. In the case of ion-momentum-recoil upon 
ionization of a model He atom, the KS momentum- 
distributions were found to be drastically wrong, dis- 
playing a single maximum instead of the characteristic 
two hump structure, and with a significantly overesti- 
mated magnitude[41]. Double-ionization probabilities 
are another example where a more sophisticated observ- 
able functional is needed [42]. 

A third problem arises when the true wavefunction 
evolves far from an SSD [10, 38]. An illustration of this 
is in certain quantum control problems, where a laser 
pulse is found to take the system to a specified target 
e.g. populating an excited state. For example, a pulse 
can be found to move the Is^ ground state of Helium 
to the l.s2p excited state. As TDDFT stays in a SSD, 
with a single doubly-occupied orbital, for all times, it 
has great difficulty in describing a state that fundamen- 
tally is best described by two SSDs. It should be empha- 
sized here that in principle, TDDFT can describe this sit- 
uation, however the exchange-correlation potential and 
observable-functionals are extremely difficult to approx- 
imate. 



B. Time-Dependent Density Matrix Functional Theory 
(TDDMFT) 

To overcome some of the difficulties of TDDFT, we 
may attempt to use the one-body reduced density ma- 



trix, defined below, as the fundamental variable. This 
has the advantage of containing more information than 
the density while retaining similar concepts and system- 
size scaling as TDDFT. An advantage is that functionals 
for the momentum distribution and kinetic energy are 
given exactly in terms of the density matrix. The first- 
order spin-summed reduced density matrix is defined 
as: 

^'*(rVi, r2CT2---rjVCTAr, t)*(rCTi, r2CT2---rjVCTAr, (2) 

and can be diagonalized by the so-called natural or- 
bitals,Cj(r,i): 

p(r',r;i)=E^^W^l(i-''^)O(r,0 (3) 

here < rij{t) < 2 is the time-dependent occupation 
number of the jth natural orbital (NO), and r]{t) = 
N. We can now interpret the quantum control example 
discussed at the end of Sec. II A as an issue with the NO 
numbers. The NOs of the true system will begin with 
one NO occupation near 2 and the rest close to zero, and 
will end with two occupation numbers close to 1. The 
KS system, beginning in the ground-state SSD, has only 
a single natural orbital which is doubly occupied and 
time-evolution in the one-body KS Hamiltonian means 
that the NO occupation number always remains at 2. By 
using the density matrix of the true system as basic vari- 
able instead, we allow the NO occupations to change 
and thus it should be better than KS at capturing the cor- 
rect behavior with simpler fimctional approximations. 

The equation of motion for the density matrix is 
known as part of the BBKGY hierarchy of equations: 

ip{r',r;t) = ( — ^ + Wcxt(r ;<) + — - Ucxt(r'; t) j /9(r', r; t) 

+ J dr2/„„(r',r,r2)p2(r',r2,r,r2;t) (4) 

where f^^{r'r,r2) = l/|r — r2| — l/|r' — r2| and p2 is 
the second-order spin-summed reduced density matrix, 
conveniently decomposed as 

P2(r',r2,r,r2;i) = n(r2)p(r', r; i) - p(r', r2; i)p(r2, r; t)/2 
+Pc(r',r2,r,r2;t) . (5) 

for a closed-shell spin-saturated system. Atomic units 
{e^ = h = nie = 1) are used throughout this paper. 
Here the terms play the role of Hartree, exchange, and 
correlation respectively. Setting pc = results in the 
TDHF equations. One could look to close this expres- 
sion by expressing pc as a functional of p and this is 
the usual approach in TDDMFT. One naturally seeks an 
adiabatic approximation in which the time-dependent 
density-matrix is fed into a p2-functional bootstrapped 
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from ground-state density-matrix functional theory [43- 
46] Unfortunately neither TDHF nor the adiabatic fimc- 
tionals lead to NO occupation numbers changing in 
time [11-13]. Thus we must look further afield to try 
and fix this. 

A possible solution to this was proposed in Ref. [10] 
where pc is given by an auxiliary semiclassical dynam- 
ics calculation running parallel to the density-matrix 
evolution. That is, we time-evolve p using Eq. 4 with 
Eq. 5, where all terms are treated exactly quantum- 
mechanically except the last term pc of Eq. 5. This last 
term is treated as a driving term in this equation, its 
value obtained from the parallel semiclassical propa- 
gation of the system. At each time, the second-order 
density-matrix, the density matrix, and the density are 
calculated in the semiclassical system, which then uses 
Eq. 5 to extract a semiclassical pc. This term is then 
inserted in Eq. 4 as a driving term. In Ref. [10], it 
was argued that such an approach addresses the main 
obstacles TDDFT and adiabatic TDDMFT approaches 
for real-time dynamics: memory is naturally carried 
along by the classical trajectories, initial-state depen- 
dence is automatically accoimted for, all one-body ob- 
servables may be obtained without the need for fur- 
ther observable-functionals, and NO occupation num- 
bers are time-evolving. We can expect this method to 
most improve upon the purely semiclassical calculation 
when in the high-density limit, as the dominating ex- 
change affects will be treated exactly. 

Various candidates for the semiclassical method to 
use in this scheme are possible, with their numerical 
efficiency inversely related to their semiclassical rigor 
and accuracy [10]. Most efficient and least accurate is to 
use a quasiclassical propagation of the one-body Wigner 
fimction [27, 28], while least efficient but most rigorous, 
solving TDSE exactly to order h, is to use the Heller- 
Hermann-Kluk-Kay (HHKK) propagator [14-16, 29-31]. 
The Frozen Gaussian (EG) method [26] may be viewed 
as an approximation to HHKK, and was argued to 
be a good candidate for use in a "forward-backward" 
scheme for pc [10]. We focus in this paper therefore on 
EG dynamics, and our present goal is to investigate how 
accurate it alone is for electron dynamics. Although in 
the scheme of Ref. [10], semiclassics is used only for the 
correlation component of the dynamics, the scheme is 
likely to be most accurate when the semiclassical prop- 
agation of the whole system is reasonably accurate. The 
purpose of the studies in Section 111 is therefore to study 
how accurate EG dynamics on model systems is. 



C. Frozen Gaussian Dynamics 

Since the earliest days of quantum mechanics semi- 
classical methods have been explored, for the purpose 
of interpreting and understanding quantum mechan- 
ics via the more intuitive classical dynamics, and also 
for the purpose of approximation. These methods con- 



struct an approximate quantum propagator utilizing 
classical trajectory information alone. Although there 
are a variety of forms (e.g. [14, 15, 26, 31, 47-51, 53]), 
their essential structure is a sum over classical trajec- 
tories: Eci.traj. C'UOe*"^''*'^'"' where S'^(t) is the classi- 
cal action along the ith trajectory, and the prefactor 
Ci{t) captures fluctuations around the classical path. 
Miller [49] showed the equivalence of different semiclas- 
sical representations within stationary-phase evaluation 
of the transformations. Semiclassical formulae have 
been derived both from largely intuitive arguments (e.g. 
Heller's frozen and thawed gaussians Ref.[26, 52]) as 
well as from rigorous asymptotic analyses of the quan- 
tum propagator (see e.g. Refs.[16, 50, 51, 53]) that sat- 
isfy TDSE to order h. The latter are based on taking 
the semiclassical limit of Eeynman's path integral for 
exact quantum dynamics. This yields a propagator of 
van Vleck form [50], that requires solving a boundary 
value problem to find the classical paths eg. from x' to x 
in time t; transforming these to "initial-value" represen- 
tations where instead a sum over an initial coordinate- 
momentum phase-space is performed, makes the nu- 
merics significantly more feasible, especially for longer 
times and more degrees of freedom. 

Semiclassical dynamics captures quantum effects 
such as interference, zero-point energy, tunneling (to 
some extent), while generally scaling favorably with the 
number of degrees of freedom. As the propagator is con- 
structed from classical trajectories, intuition about the 
physical mechanisms underlying the dynamics can be 
gained. Although mostly used for nuclear dynamics 
in molecules, and condensed phases, there have been 
a handful of applications to electrons [17-23]. The rea- 
sons for the reluctance of applying semiclassics to elec- 
trons were recently discussed in Ref. [23]. One is that the 
classical dynamics of interacting electrons are typically 
mixed with regular and chaotic regions. This makes 
the semiclassical sum over trajectories increasingly dif- 
ficult to converge as time evolves. Another is the clas- 
sical autoionization problem discussed in the introduc- 
tion; we return to this in Sec. IIIC. A third problem is 
that the classical equations of motion become singular 
at the nucleus, where the potential diverges Coulombi- 
cally, requiring the need for regularization methods to 
be applied. As Ref. [17] showed however, this is ap- 
parently only a problem when systems are treated in re- 
duced dimensionality: in three-dimensions, no special 
techniques are required to deal with the Coulomb poten- 
tial except for a set of measure zero trajectories which hit 
the nucleus head-on (and these trajectories can be safely 
discarded). Ref. [23] showed that all these obstacles can 
be overcome to make semiclassical calculations of atoms 
and molecules quite practical. 

The EG propagation we are exploring can be ex- 
pressed mathematically as a simplified version of the 
HHKK propagator [14-16, 29-31] where the A^-particle 
wavefunction at time t as a function of the 3A^ coordi- 
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nates we denote r, is: 

- / ^^{LhtPt)e^'^^'{noPom (6) 

where {qt,Pt} are classical phase-space trajectories 
in 6iV-dimensional phase-space, starting from initial 
points {qo, po}, and is the initial state. In Eq. 6, (r |qp) 
denotes the coherent state: 

(l|qp) = n (— ) e-*('-^-«^)'+^P^("^-9^-)/'^ (7) 

where jj is a chosen width parameter and St is the clas- 
sical action along the trajectory {qt, Pt}. In HHKK, each 
trajectory in the integrand is weighted by a complex pre- 
factor based on the monodromy (stability) matrix. The 
pre-factor is time-consuming to compute, and scales cu- 
bically with the number of degrees of freedom, but in 
the FG approximation, it is set to unity. As a conse- 
quence, although HHKK solves TDSE exactly to order 
h, the FG propagator does not; also, although it can 
be shown that HHKK results are independent of the 
choice of width parameter jj, FG results are not. For 
all our calculations we take ■yj = 1. Neither the HHKK 
propagation nor FG are unitary; typically we find the 
norm of the FG wavefunction decreases with time, and 
so we renormalize at every time-step. (In fact for the 
typically chaotic dynamics of interacting systems, the 
HHKK pref actor, and norm, grows with time [23]). 

The phase-space integral in Eq. 6 is done by Monte 
Carlo, with the distribution of initial phase-space 
points weighted according to the initial distribution, 
|(qoPo|^'i:)p. Although Monte-Carlo methods scale as 
\/N for positive integrands, the oscillatory phase can 
make the FG propagation very difficult to converge. The 
scheme of Ref. [10] takes advantage of the "forward- 
backward" nature of the propagation of the second- 
order density-matrix (i.e in the second-order density 
matrix, there is both a ^'(t) and a ^*(t)), which leads to 
some cancellation of phase for more than two electrons. 
We also observe that the spatial-symmetry of the initial- 
state is preserved during the evolution (since the Hamil- 
tonian is for identical particles, exchanging coordinate- 
momentum pairs of two electrons does not change the 
action). As we study here two-electron singlet states, the 
wavefunction is spatially-symmetric under exchange of 
particles. 

III. RESULTS+DISCUSSION 

We study the accuracy of FG propagation for sev- 
eral different systems involving two soft-Coulomb- 
interacting electrons in a spin-singlet in one-dimension. 
Because the TDSE for two interacting electrons in one- 
dimension can be mapped onto a TDSE for a single 
electron in two dimensions, the exact problem is eas- 



ily numerically solvable, and we can compare our FG 
results to exact ones. The exact ones are obtained ei- 
ther using the approximate enforced time-reversal sym- 
metry (aerts) algorithm coded in the octopus code [63], 
or our own exponential mid-point rule code (both use a 
fourth order Taylor expansion for the exponential). We 
also comment on how these compare with results us- 
ing various common approximations in TDDFT for ei- 
ther the exchange-correlation functional (Sec. Ill A), or 
for observable-functionals (Sees. Ill B and III D). 

In the latter case, we compute observables from usual 
operators acting directly on the exact KS wavefunc- 
tion: this enables us to isolate errors arising from the 
observable-functional approximation alone (i.e. as op- 
posed to errors from the approximation used for the 
exchange-correlation potential). For two electrons in 
a spin-singlet, given the exact time-dependent density, 
it is straightforward to obtain the exact KS wavefunc- 
tion [36, 54-56]. We briefly review how. In one- 
dimension, the continuity equation reads 

dn{x,t) dj{x,t) 

' a ~ ^ ' 

at ox 

and since the exact KS system reproduces the exact den- 
sity at all times, 

j{x,t) =jKs(x,t) = 2$ne(/)Ks(a;,<)- ^ • (9) 

The last equality follows from the fact that for two elec- 
trons in a spin-singlet, there is just one spatial KS orbital 
0KS (doubly-occupied). Writing ^ks as an amplitude 
times a phase, we deduce, 

t) = /^exp (^ f ^c/x') . (10) 

So, the procedure is as follows: we solve the TDSE ex- 
actly numerically, finding the exact two-electron corre- 
lated wavefunction. The density and current are then 
extracted from n{x,t) = 2\(l)Ks{x,t)\'^ and Eq. (8), and 
then inserted into Eqs. (10). Usual approximations 
for observables in TDDFT evaluate the operator corre- 
sponding to the observable directly on the KS state; al- 
though this is exact for observables directly related to 
the density, there are errors for other observables, e.g. 
the momentum-density, which is what we shall look at 
in examples III B and III D. 

In all cases, we will utilize the exact ground-state 
wavefunction to construct our initial state for propaga- 
tion. This is to separate out any errors due to a poor 
ground-state. For a general case where the initial wave- 
fimction is unknown exactly, an appropriate approxi- 
mate initial state (e.g. sum of a few KS SSD's, or from 
a high-level wavefunction calculation) would be used. 
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A. Hooke Model: Spectra 



We begin with the Hooke model for a one- 
dimensional two-electron quantum dot defined by the 
Hamiltonian 

The Hamiltonian becomes separable when written in 
relative r — x\— X2 and center of mass R = [xi + X2)/'2. 
coordinates, and we use the quantum numbers of these 
systems to label the excitations of the full system. In 
Table I, we give the frequencies of the lowest 8 singlet 
excitations of this system, along with labels in i?, r. The 
r label must be even in order to have the correct spatial 
symmetry for a singlet state, while even/odd values for 
R determine the parity of the state. We have grouped 
the levels into multiplets whose members become de- 
generate if the electron-interaction is turned off. Each 
of the two states in the second and third multiplets are 
states of double-excitation character: they are mixtures 
of a single excitation (where, in a non-interacting ref- 
erence, only one electron is promoted from its ground- 
state orbital to an excited one), and a double-excitation, 
where both the electrons occupy excited orbitals [34, 35]. 
The three states in the fourth multiplet are mixtures of a 
single-excitation and two double-excitations. 

Such states of double-excitation character are well- 
recognized to be poorly described with the adiabatic 
approximation in TDDFT, which is used in almost all 
calculations [33, 34]: instead of producing a multiplet 
of states, only one state is obtained in the adiabatic ap- 
proximation, and it has only single-excitation character. 
In order to find the correct frequencies, TDDFT requires 
a frequency-dependent kernel to mix the KS single and 
double-excitations, as was derived in Ref . [34] . The ker- 
nel is an a posteriori correction to an adiabatic approxi- 
mation and although it has been recently tested on a va- 
riety of systems [57-60], it is not widely used; almost all 
codes still rely on the adiabatic approximation that fails 
to capture these. This is the motivation behind our first 
exploration of the FG dynamics: is semiclassical corre- 
lation able to capture states of mixed single and double- 
excitation character? The answer, as shall see shortly, is 
yes, at least approximately. 

We compute the FG spectra by Fourier-transforming a 
time-propagation. It is common to start in a linear kicked 
state defined as 

«'fc(xi,X2) = e'''^^'+^''>^o{xi,X2) (12) 

where is the exact ground-state wavefxmction, in or- 
der to populate the excited states and then to Fourier 
transform the time-dependent dipole moment, d{t), to 



TABLE I: The singlet excitation frequencies ujn = En — Eo, 
where the ground-state energy Eo — 1.774040 a.u., for the 
Hooke model. For each excitation, the energy can be written 
En = Er + Er, where R,r are labels in the center of mass and 
relative coordinates respectively. 



Label (R,r) 




, ,FG 


1,0 


1.000000 


1.0 


0,2 


1.734522 


1.6 


2,0 


2.000000 


2.0 


1,2 


2.734522 


2.6 


3,0 


3.000000 


3.0 


0,4 


3.648334 


3.5 


2,2 


3.734522 


3.6 


4,0 


4.000000 


4.0 



find the dipole-spectra: 

d{uj) = J dte'"*d(t) (13) 

The kick is equivalent to applying an impulse electric 
field, Svcxtix, t) = kS{t)x, to the system[61-64]. For The 
spectra found this way consists of peaks at the excitation 
frequencies, the widths of which becomes smaller the 
longer you propagate the system. For Hooke's model, 
this method needs a slight adjustment. In linear re- 
sponse, the kick yields a perturbation that is propor- 
tional to R, which only couples the ground state (0, 0) 
to the (1,0) state, as the Hamiltonian is simply a har- 
monic oscillator in this coordinate; the harmonic poten- 
tial has the special property that a dipole can only con- 
nect two states differing by one quantum number. This 
can also be seen as a consequence of the Harmonic Po- 
tential Theorem[65] . So in the exact system, we only see 
a single peak in the dipole moment when using a linear 
kick. In order to see higher excitations, we use quadratic 
and cubic kicks and plot the quadrupole and third-order 
moments respectively. For example the quadratic kick, 
gjfc(£c-+y-)^ in the linear-response regime will be propor- 
tional to both R^ and r^, which, using the symmetry 
and parity discussion earlier, will couple most strongly 
to the (2, 0) and (0, 2) states. 

We begin the FG calculation in the quadraticaUy kicked 
state (using the exact ground state) with k = 0.01 and 
use 120000 classical trajectories. These trajectories are 
propagated using the standard leapfrog verlet algorithm 
for a total propagation time of T = 200 a.u with the 
dipole moment calculated every 0.1 a.u. The power 
spectrum (absolute value squared) of the quadrupole 
moment is shown as the solid line in Fig. 1; we use the 
power spectrum as it is usually more stable than either 
the purely real or imaginary part, and we also apply the 
3^'' order polynomial filter given in Ref. [62] to reduce 
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noise when computing the Fourier transform. We ob- 
serve two peaks, representing each of the two states in 
the second multiplet ((0,2), (2,0)). The FG therefore does 
capture states of double-excitation character. The exact 
values of the (0,2) and (2,0) excitations are indicated on 
the graph as dashed lines. We see that FG gives the peak 
at cj = 2 a.u. exactly right and the other is shifted lower 
to approximately lo = 1.58 a.u.. The peak at 2 a.u. is 
purely an excitation in the center of mass coordinate, 
where, as mentioned, the system is purely harmonic. 
For harmonic systems, most semiclassical methods, in- 
cluding FG, are exact [26, 50], and so it is expected that 
FG works well for this peak. 

We reiterate that TDDFT in its usual adiabatic approx- 
imation, only yields one peak; Ref. [35], examined this 
pair of excitations using adiabatic exact exchange within 
the single-pole approximation and found, as expected, a 
single peak (of frequency 1.87 a.u.). 

A quadratic perturbation yields zero transition prob- 
ability to states of odd parity, such as in the third mul- 
tiplet. To access these states, we cubicly kick the system 
and find the third-order moment, also plotted in Fig. 1 
(dashed line). This is a harder test for FG as we are am- 
plifying any errors far from the center in the decaying 
part of the wavefunction. Again we see that FG works 
well, giving three peaks corresponding to the first exci- 
tation, and the mixed single- and double-excitations in 
the third multiplet, with again the center-of-mass exci- 
tations given exactly. Looking at Table I, we see that 
the FG correctly reproduces the property that the pair 
of peaks in the third multiplet differs from those in the 
second by a single excitation in the center-of-mass coor- 
dinate. 

The peaks in the fourth multiplet are most clearly re- 
solved by applying a quartic kick to the system. These 
three states involve mixtures of a single-excitation with 
two double-excitations, and FG provides a good approx- 
imation. 

Our spectra could be made cleaner by running for 
longer times or using harmonic inversion methods for 
example. However they are adequate for our current 
purposes of illustrating that the FG method captures 
correlated states of mixed single and double-excitation 
character. 

In conclusion, although the original motivation for a 
semiclassical description of electron correlation was to 
address challenges TDDFT has for real-time dynamics 
in non-perturbative fields [10], rather than for spectra, 
this study shows that semiclassical dynamics may nev- 
ertheless also be useful in the linear response regime. In 
particular, FG dynamics does capture states of double- 
excitation character, albeit approximately, missing in the 
usual adiabatic approximation of TDDFT. In compar- 
ison to the "dressed TDDFT" that uses a frequency- 
dependent kernel derived in Ref. [34], the FG results 
are not as accurate, e.g. for the ((0,2),(2,0)) multiplet, 
dressed TDDFT gives 1.712 a.u. and 2.000 a.u. On the 
other hand, in the dressed TDDFT approach, the proce- 
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FIG. 1: The power spectra of the Frozen Gaussian quadratic, 
cubic, and quartic moments for the Hooke's model dot. The 
positions of the exact frequencies are shown as the vertical 
dashed lines. 



dure to obtain the double-excitations involves an identi- 
fication by the user of zeroth-order single- and double- 
excitations that are likely to interact strongly, while in 
the present semiclassical approach, they emerge natu- 
rally from the dynamics. 

We make two remarks at this point. First, a full HHKK 
treatment would likely yield more accurate semiclassi- 
cal spectra, given that the prefactor missing in FG is 
complex, so its phase contributes interference effects im- 
portant in determining the resonant frequencies. Sec- 
ond it is important to note that the proposed method 
of Ref. [10] treats only the correlation component of the 
density-matrix dynamics semiclassically, rather than of 
the entire djmamics as we have done here; this may re- 
sult in improved accuracy of the spectra and we shall 
investigate this in the future. Given the more favorable 
scaling with system-size of FG than of HHKK, and given 
that it will ultimately be used in conjunction with exact 
kinetic, Hartree, and exchange components of the evo- 
lution, it is the FG method we are most interested in at 
present. 



B. Hooke Model: Natural Orbitals and Momentum 
Distributions 

Next, we apply a resonant driving perturbation to 
the Hooke model dot, to investigate whether FG yields 
time-dependent NO occupation numbers accurately. 
Due to the harmonic potential theorem [65], applying 
an electric field does not change the natural orbital oc- 
cupation numbers; instead we apply a quadratic per- 
turbation, 5vcxt{x, t) = k{t)x'^ /2, with a time-dependent 
spring constant: 

fc(t) =0.05sin(2i) (14) 
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tribution, defined as 



FIG. 2: NO occupation numbers of the Hooke model dot 
driven at a resonant frequency of 2.0a.u. 




FIG. 3: The quadrupole moment for Hooke model dot driven 
at resonant frequency. Due to its highly oscillatory behavior, 
we show the exact values only as square points for clarity. 



designed to encourage population transfer to the (2, 0) 
state. 

In Ref. [10], the NO occupation numbers for 
the Moshinsky model were calculated, however the 
Moshinsky Hamiltonian is completely quadratic and so, 
as mentioned, FG is exact. In the Hooke model case, the 
soft-Coulomb interaction between the electrons makes 
it a little more realistic model of a quantum dot, and a 
harder test for FG. In the following FG calculations for 
the system being driven with Eq. (14) , only 60000 clas- 
sical trajectories were needed for convergence. 

In Fig. 2, we plot the exact and FG NO occupation 
numbers. Recall that for both exact TDDFT and all com- 
mon approximations in TDDMFT, the occupation num- 
bers are fixed by their initial values and cannot change. 
It is quite clear that despite not being perfect, FG works 
very well here, tracking the large changes in the occupa- 
tion values. The quadrupole moment is shown in Fig. 3, 
showing that the FG density is also very accurate. 

In an exact TDDFT KS calculation, the density would 
be exact, but the NO occupation numbers would remain 
at integers 2 (occupied orbital), and for all the rest. We 
now show how this feature drastically affects observ- 
ables that are not directly related to the coordinate-space 
density. We shall focus on the one-body momentum dis- 



n{p;t) 



dp' 



d.Tdx'e*(P^+f'^')4'(a;',a;;<) 



(15) 

i.e. it is the probability of finding any one electron with 
momentum p. Although it is known that the momen- 
tum distribution of the exact KS wavefunction is not that 
of the true system [9, 41], little is known about how to 
extract the latter from the former. In the absence of a 
good observable-functional for momentum, one resorts 
to simply using the exact KS momentum distribution. 

In Fig. 4, we show the exact, exact KS (computed di- 
rectly from the exact KS orbital obtained via the proce- 
dure in Sec III), and FG momentum distributions at four 
snapshots in time. The first moment of the KS and ex- 
act distributions must integrate to the same value for all 
one-dimensional systems (in this case zero), due to the 
facts that (p(t)) — J dxj{x,t) = J dxxh{x,t), and that 
the KS exactly tracks the density. Despite this, the KS 
distribution oscillates wildly. The reason for this is not 
dissimilar to that seen in Ref. [9]: the KS system becomes 
strongly non-classical as it attempts to describe the sys- 
tem of two-electrons using a single (doubly-occupied) 
orbital. For short times, the KS n{p) behaves well, how- 
ever as the breathing motion of the system becomes 
more pronounced, moving out further and faster, it be- 
comes an increasingly non-classical dynamics for a sin- 
gle orbital, such as in KS, to describe. The underly- 
ing phase-space distribution of the KS system develops 
strong oscillations into negative regions, signifying the 
non-classical description of different "parts" of one elec- 
tron moving in different directions. The momentum- 
distribution represents one observable that the occupa- 
tion numbers strongly influence. The FG results are 
much more accurate than the exact KS, as is also re- 
flected in the kinetic energies shown in Fig. 5. 



C. Soft-Coulomb Helium: The Problem of Classical 
Autoionization 

Now we move from model quantum dots to model 
atoms: a soft-Coulomb interaction is used for the nu- 
clear potential, as a one-dimensional model of a Helium 
atom (see e.g. Refs. [56, 66-69]). The Hamiltonian is 
given by: 



H 




^/{xl-X2)^ + l 



(16) 



where the third term on the first line represents an ex- 
ternal electric field applied to the system. This model 
has been used in several studies of correlated electron 
dynamics in strong fields as well as in TDDFT, since it 
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FIG. 4: The exact (solid line), KS (short-dashed line), and FG 
(long-dashed line) momentum densities calculated at times 75, 
135, 160, and 200 au for the on-resonant hooke system. 




FIG. 5: The kinetic energy for the exact, KS, and FG calcula- 
tions, when the Hooke dot is driven on resonance. Although 
the KS system must reproduce the exact momentum, the ki- 
netic energy becomes worse for longer times, while the error 
in the FG remains about the same. 



allows a numerically exact solution to be calculated rel- 
atively easily. 

With our electronic system no longer everywhere 
bound, we encounter the problem of classical autoion- 
ization, mentioned in the introduction. Even when no 
field is present, this haunts the dynamics, as we will 
now show. We begin in the exact ground-state of the 
Hamiltonian Eq. (16) with e{t) = 0. Fig. 6 plots the 
positions xi,X2 of the trajectories in the initial distri- 
bution, and where they have evolved to at T = 10a. u. 
The initial distribution centered at the origin increas- 
ingly evolves into a cross, signifying the classical au- 
toionization: one electron zipping off to infinity after 
stealing the energy from the other left near the origin. 
Using just 20000 trajectories to compute the FG sum, we 
plot in Fig. 7 the density at four snapshots in time. Au- 
toionizing wavepackets are clearly seen, evolving out 
quasiclassically on either side of the central distribu- 



tion. (Note that in the exact problem, nothing happens; 
the distribution remains centered at the origin, as the 
initial state is an eigenstate.) Although it may appear 
that the autoionizing packets are growing as a function 
of time, this is an artifact of the way the FG dynamics 
is renormalized (see Sec. II C). Before renormalization, 
the FG norm decreases with time as more trajectories 
from the central distribution are lost to autoionization 
and other incoherent effects. Why some of the classi- 
cally autoionizing trajectories seem to add coherently, 
forming the wavepackets on the sides, remains to be 
understood. However choosing a different initial seed 
to generate the initial distribution results in autoioniz- 
ing wavepackets centered at different positions moving 
with different speeds. That is, certainly the results in 
Fig. 7 are imconverged. In fact, in Fig. 8 we show the ef- 
fect of increasing the number of trajectories is to signifi- 
cantly decrease the classical autoionization. We deduce 
that by phase interference in the semiclassical sum, con- 
tributions from classically autoionizing trajectories can- 
cel each other out. Considering that HHKK-semiclassics 
gives the correct dynamics to order h, this is as it should 
be: Although FG dynamics is not correct to order h a 
similar effect happens here. 

Although these results show that classical autoioniza- 
tion is not a fundamental problem of the semiclassical 
method, but rather a convergence issue, it does remain 
an important practical one for most applications; impos- 
sibly large numbers of trajectories could be required to 
make the classical autoionization effect small enough 
for the purposes. In our following studies we termi- 
nate trajectories in one of two ways. In the applications 
we consider, we do not expect significant amplitude be- 
yond a certain distance. So, (i) we simply terminate any 
classical trajectory where the coordinate of one electron 
has reached a prescribed distance, L, away from the nu- 
cleus. This is similar to what is done in Ref. [23] in ob- 
taining semiclassical spectra (where the termination dis- 
tance was additionally dependent on the total energy of 
the trajectory). Another alternative is to use purely the 
energy terms to determine autoionization such as was 
done in Ref. [70]. At all times, this means that we may 
still observe autoionizing wavepackets peeling off from 
the central distribution, that eventually reach L, only 
after which they are discarded; therefore at the earlier 
times, these may erroneously contribute to the FG. We 
therefore also consider (ii) prescreening all trajectories, 
such that if they, at some point during the duration of 
interest, reach L, they are discarded from the very start 
of the propagation. (Note this prescreening process is 
very fast, as only classical d5mamics needs to be run; 
therefore many initial candidate trajectories may be ex- 
plored). In this way no autoionizing trajectories (as de- 
fined by L) ever contribute to the FG sum. In our partic- 
ular studies, there was some difference between the re- 
sults obtained with (i) and (ii) especially at shorter times, 
with (ii) leading to narrower distributions than (i). It is 
imclear to us which way is more "correct", especially 
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FIG. 8: The FG density at time t — 25 a.u. (with no field ap- 
plied to the soft-Coulomb He atom). Adding more trajectories 
reduces the classical autoionization. 



FIG. 6: The positions, xi and X2, of both particles plotted 
against eachother at the initial time (square points) and at time 
t = 10 a.u. (triangle points). The problem of classical autoion- 
ization is characterized by the cross-shape distribution at the 
later time: one particle falls to the center while the other flies 
away from the atom. 
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FIG. 7: Snapshots of the FG density at various times when no 
field is applied to the soft-Coulomb He atom. The lobes on 
either side of the central distribution represent autoionizing 
wavepackets. The number of trajectories in the monte carlo 
sum was M = 20000. 



given the fact that, theoretically, none of the trajectories 
are "wrong": including all, together with their phases, 
there is no classical autoionization problem. In particu- 
lar, does the earlier behavior of a trajectory that eventu- 
ally autoionizes, contribute in a non-cancelling, sensible 
way to the dynamics at times prior to its autoionizing 
event? This is an important question we leave for future 
work. 



D. Soft-Coulomb Helium:Strong field 



We now apply a relatively strong field to the model 
He atom: 



e(t) = ^ cos(0.5i) 
v2 



^ t<20 
1 t>20 



(17) 



where the electric field is ramped-up linearly until t ~ 
20 au. The value of this field is large enough that classi- 
cal autoionization effects are relatively small in compar- 
ison to the large oscillations in the dynamics induced by 
the field. We use the procedure (i) described above and 
terminate trajectories once |a-i| > 30a.u. We found the 
results were essentially converged starting with 2 x 10^ 
trajectories and finishing with 5 x 10^, the rest being dis- 
carded by the termination condition at some point dur- 
ing the run. We plot the NO occupation numbers in 
Fig. 9, up to a time of 35 a.u. Notice that the number of 
trajectories needed for convergence in the soft-Coulomb 
He atom is much larger than those needed in the Hooke- 
model quantum dot; the latter is almost a best case sce- 
nario for semiclassics, because of the quadratic nature 
of the external potential. The FG captures the change 
in occupation numbers, perhaps a little too enthusias- 
tically; Figure 10 shows the dipole is also reasonably 
well approximated. In Fig. 11, we plot the momen- 
tum distributions at several snapshots of time, compar- 
ing the FG with the exact, and with the exact-KS, as in 
Sec. IIIB. We see that the error of the FG remains about 
the same throughout the evolution, while the error of 
the KS increases with time. This is a direct consequence 
of the change in NO occupation numbers: the KS mo- 
mentum distribution is that of a SSD with NO occu- 
pation numbers 2 for one orbital and zero for all oth- 
ers, while that of the true and, roughly captured by FG, 
changes dramatically in time, as shown in Fig. 9. Con- 
sequently, typical observables not directly related to the 
density, when evaluated using the usual operators on 
the KS wavefunction become badly approximated; the 
imknown observable-functionals of exact TDDFT are 
likely quite complicated. 
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FIG. 9: NO occupation numbers for the soft-Coulomb He atom 
dynamics induced by Eq. 17. 
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FIG. 10: Comparing the FG dipole moment in the soft- 
Colulomb He dynamics under Eq. 17. 



E. Soft-Coulomb Helium: Towards Optimal Control 
Theory 

Our final example is also the hardest one: we ap- 
ply an optimal control field to the soft-Coulomb He 
atom, to try to get it to evolve from the ground-state to 
the first-excited state. The occupation numbers in the 
ground-state are close to 2 for the highest, and close to 
zero for all others, while in the target state, the high- 
est two occupation numbers are both close to 1, having 
a double Slater determinant character. (See also Sec. II, 
and [10, 38]). However, instead of optimizing the field 
for the FG evolution, we simply take an optimal field 
found for the exact dynamics, and ask how well does it 
work for the FG problem. (Of course, the former is what 
we ideally would do, however, we leave the coding of 
optimal control problems for future work.) We find the 
optimal field using the optimal control functionality of 
the octopus code [63, 71]. Moreover, we consider a rela- 
tively short duration for the optimal pulse, T = 35 a.u., 
as we will find this illustrates already the challenges the 
FG approach has for such problems. Because there are 
only a few cycles of the laser field in this short time, the 
yield in the exact problem is not very large: the optimal 
field found in Fig. 12 yields a final state population of 
0.73.. The field is significantly weaker than in the pre- 
vious section, and we really do not expect much ioniza- 
tion, classical or real. To avoid the situation of Fig. 7, we 
therefore use method (ii) to terminate the trajectories. 




FIG. 11: The momentum density n(p) at times 10, 15, 20, 25, 30 
for the strong field case. Solid line is exact, dotted is KS, and 
short-dashed is FG. 
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FIG. 12: The optimal electric field to maximally populate the 
first excited singlet-state of soft-coulomb Heliimi within a time 
of 35au. 



prescreening them such that no trajectories where one 
coordinate reaches a distance L = 10 a.u. at any time in 
the run are included. We find that beginning with 5x10^ 
candidate trajectories, the prescreening process discards 
many such that 2391950 remain for the FG calculation. 

The result for the NO in in Fig. 13 are disappoint- 
ing. The FG results are converged, as different seeds, 
and increasing the number of candidate trajectories did 
not change the results very much. Also, simply rimning 
unscreened calculations with increasing numbers of tra- 
jectories, appeared to converge, more or less, to these re- 
sults. Although initially the largest FG NO occupation 
number tracks the exact reasonably well, we see that af- 
ter about T = 15, it turns upwards instead of continuing 
down; the FG state appears to reduce to a SSD instead 
of developing more of a double -SD character. Turning 
to Figs. 14 and 15, we see that after beginning to spread, 
the density then contracts and settles into a narrow dis- 
tribution in the well. The dipole at shorter times is not 
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FIG. 13: The highest natural orbital occupation for a pre- 
screened FG calculation with the optimal electric field of Fig. 
12. The exact NO value is nearly 1 at the end, reflecting the 
double ssd character of the target state. 



bad but fails at larger times. We also note that the den- 
sity appears to tighten as time evolves; this is a conse- 
quence of autoionization and occurs even if we do not 
apply any termination to the trajectories, in the limit of 
convergence; the predominant trajectories that autoion- 
ize are those in the tails of the distribution. 

Why the FG results are poor at larger times is, we be- 
lieve, due to the FG resonant frequencies being so off-set 
from the true ones, that a field that is optimal for popu- 
lation transfer for one is off -resonant for the other, lead- 
ing to little transfer. The exact excitation frequency to 
the first excited state is 0.53, while the FG one is about 
0.68 (found by applying a kick to the initial ground- 
state, as in Sec. Ill A). For example, in the exact case, 
if the frequency of the applied field is shifted by as little 
as O.OSau away from resonance, there is no longer the 
strong change in NO numbers that we saw in Fig. 2. 

Thus the optimal control problem presents an ex- 
tremely tough test of the frozen gaussian dynamics as 
it depends critically on subtle interference effects of on- 
resonance oscillations. As stated in Sec. Ill A, we expect 
the FG excitation frequencies to improve when coupled 
with the TDDM propagation, and hence also improve 
the optimal control results. This will be investigated in 
future work. 



IV. CONCLUSIONS AND OUTLOOK 

The semiclassical FG method provides an intuitive 
picture of quantum mechanics, as based on classical tra- 
jectories, each smeared out into a little fuzzy ball in 
phase-space, evolving classically in time, and added 
together with phases determined by their classical ac- 
tion. It, and its more rigorous HHKK version, has also 
provided useful numerical tools in quantum molecu- 
lar d5mamics. Its application to interacting electronic 
systems has been hesitant, although recent work [23] 
clearly shows its promise for electronic spectra. The re- 
sults of the present work suggest that for systems of in- 
teracting electrons in external fields, FG dynamics can 




FIG. 14: The density of the prescreened FG calculation with 
the optimal field at various times compared to the exact. 
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FIG. 15: The prescreened FG dipole moment versus the exact 
as a function of time. Although initially the FG mimics the true 
at short times, the delicate nature of quantum control means 
FG is too inaccurate for this problem. 



also be useful. 

In relation to TDDFT and adiabatic TDDMFT, we 
have shown that FG dynamics overcomes some of 
the problems these methods have; capturing states of 
double-excitation character, changing occupation num- 
bers, and accurate momentum distributions. There are 
nevertheless issues related to its convergence, efficiency 
and accuracy, that require further study. The issues 
mainly concern the computational details of the classi- 
cal trajectories, in particular, the large number of trajec- 
tories needed for convergence and also the issue of clas- 
sical autoionization. 

Although a large number of trajectories is required 
for convergence for our two-electron model atoms, we 
expect that the "forward-backward" nature of the semi- 
classical computation in Ref. [10] for the two-body re- 
duced density-matrix will lead to favorable scaling with 
the number of electrons. Because the prescription of 
Ref. [10] only calls for the correlation component of the 
dynamics to be treated semiclassically, less accurate and 
more efficient semiclassical methods, such as thawed 
Gaussian or even just quasiclassical propagation (evolv- 
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ing an initial distribution just classically) may be worth 
exploring: the aim is just to capture "enough" correla- 
tion in a physical way to be used in conjunction with the 
exact treatment of the one-body terms in the equation of 
motion for p{r, r' ,t). 

How to deal with classical autoionization in a con- 
sistent and practical way needs further study. When 
a strong field is present and electrons are unbound at 
least for some time, the effect of classical autoionization 
may be relatively small enough for times of interest that 
a simple termination of trajectories when they get "too 
far away" makes sense. However in other situations 
prescreening trajectories by discarding at the start those 
which at any time reach a given boundary, may be a bet- 
ter approach to avoid autoionizing trajectories "on their 
way out" from distorting earlier dynamics. Yet there is 
a delicate balance: if prescreening is done over too long 
time, too many trajectories that are important at smaller 
times get discarded. The solution of this problem also 
depends on what is the quantity of interest, e.g. whether 



it is time-resolved densities, or time-averaged spectra. 

By including correlation semiclassically, the scheme 
of Ref. [10] is likely to provide more accurate dynam- 
ics than the FG treatment of the entire dynamics shown 
here, as well as improving over the TDHF and adiabatic 
TDDMFT results. The present results are encouraging 
for this next step. 
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